Precursory dynamics in threshold systems 
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A precursory dynamics, motivated by the analysis of recent experiments on solid-on-solid friction, 
is introduced in a continuous cellular automaton that mimics the physics of earthquake source 
processes. The resulting system of equations for the interevent cycle can be decoupled and yields 
an analytical solution in the mean-field limit, exhibiting a smoothing effect of the dynamics on the 
stress field. Simulation results show the resulting departure from scaling at the large-event end of the 
frequency distribution, and support claims that the field leakage may parameterize the superposition 
of scaling and characteristic regimes observed in real earthquake faults. 
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There are many examples in nature of systems whose 
internal dynamics have as an essential feature the sep- 
aration of time scales into a long-term loading process 
and a short-term discharge. In the former, some exter- 
nal source slowly drives the variables that control the 
dynamics until these fields reach some threshold value. 
When this happens, the fast dynamics of the discharge 
takes over and generates a sequence of internal rearrange- 
ments, called an avalanche in this context, which eventu- 
ally removes the excess through internal dissipation, or 
by flow across the system's borders. These dynamical 
threshold systems have been the object of much recent 
interest, due mainly to the fact that they can be recast 
into a discrete time-evolution form. The resulting dis- 
crete dynamics transforms the fast time scale discharge 
cycle into a sequence of well-defined time steps, in each 
of which the dynamical fields undergo a stochastic trans- 
formation that depends on their configuration in the pre- 
vious step as well as annealed noise. As such, and since 
the fields are in general real valued, we will call these 
representations stochastic continuous cellular automata 
(SCCA). 

Physical systems that have been modeled and studied 
along these lines are abundant in the literature, ranging 
from the firing behavior of neural networks jjj to the 
dynamics of domain walls in magnets HI , the motion of 
vortex lines in type II superconductors |3|] , the depinning 
transition in the growth of interfaces in random media 
H , the dynamics of the energy release in solar flares || , 
and to the source processes responsible for earthquakes. 

The existence of an interevent dynamical cycle in the 
latter is suggested by recent laboratory experiments ad- 
dressing issues of solid-on-solid friction. A stable slip, 
with a slow velocity that increases with the stress level, 
is observed prior to failure, leading to a partial release of 
the accumulated stress This stress leakage mecha- 

nism is analogous to a temperature-dependent viscosity 
that has been observed in laboratory for the creeping of 



crystalline rocks, and can be modeled by the equation 
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where s(i) and a(t) are the displacement and the stress 
at time t, an is some residual stress value to which the 
system decays after <j(t) reaches a failure threshold erp, 
K is the elastic stiffness and a measures the intensity of 
this stress leakage effect. This parameter, with dimen- 
sions of inverse time, will in general be dependent on 
both the stress level and the temperature. An analogous 
leakage mechanism has also been suggested in the con- 
text of integrate-and-fire neural networks 0], and it is 
likely that the results reported here will also hold in that 
context. 

As a summary, we will show that the introduction of 
a stress leakage process as an interevent dynamics in a 
SCCA model for earthquake faults changes in a dramatic 
way the space-time patterns it generates. In particular, 
the a value for a fault may determine its overall behav- 
ior as of a scale-invariant type or a nucleation type, with 
a mixed composition in between, reproducing features 
observed in real faults j8|. The importance of this new 
parameter in earthquake source models has in fact been 
recently evaluated. Its tuning to match the characteris- 
tics of each segment in a complex computer representa- 
tion of the fault network of southern California allowed 
the generation of space-time patterns of rupture of un- 
precedented realism || . 

We will consider here the Rundle- Jackson-Brown 
(RJB) SCCA model for an earthquake fault (To), in its 
uniform long-range interaction, mean-field version 
Extensive work has recently focused a near-mean-field 
version of this model, where the interaction range has a 
cutoff and the model can be mapped onto an Ising-like 
Langevin equation Its dynamical variables are two 
continuous real- valued fields, slip Si(t) and stress <Ji{t), 
defined on the sites i of a lattice. A constitutive equation 
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couples these fields, 



<n(t) =J2T ijSj (t)+K L vt 
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where Ty is the interaction matrix, or stress Green's 
function, Tu = -K, K = K L + K c , K c = Ej^Tij, 
Kl parametrizes the loading interaction, and v fixes a 
relation between the short and long time scales of the 
physics of the model. Equation (||) specifies the stress 
at each site as long as max[<7,(£)] < <jp, which defines 
the stress threshold o>, taken as uniform over the lat- 
tice. The interevent, or loading, dynamics of the stan- 
dard RJB model is thus very simple; the slip field remains 
static while the stress field undergoes linear growth. As 
soon as the stress at one site reaches the threshold, a fast 
stochastic relaxation dynamics takes over. The slip at 
a site that fails is discontinuously reset to a new value 
that leads the stress to assume its residual value, usually 
with some noise, introduced to represent the disorder in 
the rheology. The stress drop is redistributed among the 
interacting neighbors, with an intrinsic dissipation mea- 
sured by the factor S = K^jK: 



Act, = -£|Ao-i 
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where i is the site that failed. This increase may cause 
other sites to fail as well, and the process continues until 
all sites have stress below failure. This cascade of failures, 
or avalanche, is the model's equivalent for an earthquake. 

The relaxation dynamics of the model can be cast into 
a single field formulation that is specially convenient for 
computer simulations. In units of the short time scale, 

*i(t + 1) = «r,(t) + {Tu * 3 '®'"* *) fo(*) - °f) 
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where &(t) is the noise term. We direct the reader to 
Ref. for a complete and pedagogical discussion of 
this relaxation dynamics. 

We report in this paper results obtained for the intro- 
duction of stress leakage, as modeled by Eq. ([!]), as an 
interevent dynamics of our SCCA model. The equations 
for these intervent dynamics are obtained by taking the 
time derivative of Eq. (0), substituting dsi(t)/dt from 
Eq. ([!]), and choosing, with no loss of generality, ctr = 0, 
to get, 
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a set of N coupled equations for the stress field. Without 
introducing any particular form for the Ty, these equa- 
tions can be combined to derive the time evolution of the 
average stress for short times a(t) = <Ji{t)/N : 



da(t) 
dt 



-a-j£-a(t) + K L v, 
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with the solution, cast in dimcnsionlcss form by defining 
t = of I Klv, r\ = a I of, t — t/t, and <p = at, 
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which is an increasing (decreasing) function of dimen- 
sionless time if fj(0) < (>)1/S(f> 

Equations |^ and || above define the generic RJB model 
with stress leakage for the interevent dynamics. The cur- 
rent focus of research is on the long-range versions, where 
the interaction matrix is nonzero for a substantial 
fraction of all pairs (ij). Ideal elasticity indicates that 
this matrix should involve a dependence on the distance 
as 1/Y 3 , as expected from the results of simple labora- 
tory experiments. Nevertheless, it has been shown that 
the upper critical dimension for this interaction is d u = 2 
|]l3f and we can recover the long wavelength physics of 
this model by a simpler mean-field formulation. We will 
thus focus on the mean-field RJB model, with a uniform 
interaction matrix Ty = K c / (N — 1), where TV = 1/A is 
the number of sites in the lattice. In this case, the set 
of equations in (J5J) can be decoupled and yield a solution 
which reads, in dimensionless variables, 
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we can see that this time evolution of the stress field 
is order-preserving, i.e. r?i(0) > rjj(0) =>■ rji(t) > rjj(t), 
allowing for the determination of the initiator of the next 
event by finding the site with maximum stress at the 
completion of the previous event. The important effect 
of this leakage stress dynamics in the pattern of failures 
of the system comes from the reduction it causes in the 
statistical spread of the stress field with time. A simple 
measure of this smoothing is obtained through the time 
evolution of the variance of the stress field 



var [rj(t)] = var [rj(0)] exp 
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Because of this exponential smoothing of the stress field, 
the probability of a site to reach failure after receiving 
a transfer from a failing neighbor increases, and is an 
increasing function of the time to failure. As a conse- 
quence, the branching ratio, defined as the average num- 
ber of failures caused by each failing site, also increases. 



2 



The system is more likely to undergo larger avalanches, 
which may even be system wide when the time to failure 
is large enough. 

Because this interevent dynamics preserves order, as 
already mentioned, its introduction in a SCCA is rather 
straightforward. The site with stress nearest to fail- 
ure r\ m ax after an event is the initiator of the next, 
and solving for r in Eq. (||) for this site determines 
the time-to- failure tf, that will be used again to up- 
date the stress field over the lattice prior to the event. 
To avoid the time-consuming solving of transcendental 
equations at each time step in the simulation, the time- 
evolution equations may be linearized, as long as max 
{54>t : [(1 - <5A)/(1 - A)}4>t} « 1, which is true for our 
parameters, to read 
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and the solution for the time-to-failure is 
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FIG. 1. The frequency distribution is shown for dimension- 
less leakage factors <j> = 1.0 and 2.0; the plot with no leakage 
is shown for a comparison. The simulations were performed 
on a 128 x 128 lattice, with a dissipation factor 5 = 0.01 and 
a noise amplitude of 0.5. Data for this and all subsequent 
plots was collected from 3000000 events, after a transient of 
the same order, and was logarithmically binned. 

Figure (Q) shows a log-log plot of the frequency dis- 
tribution of events as a function of their size. The ef- 
fect of the stress leakage dynamics shows up clearly in 
the excess over scaling obtained for large events as <fi is 
increased from 0, together with a depletion of the distri- 
bution in the intermediate size range. The slope of the 
scaling part of the plot also gradually increases, starting 
from the mean-field value r = 1.5. The smoothing ef- 
fect of the leakage dynamics, together with the resulting 
larger stress average that it causes in the field as a whole, 
results in a higher probability for large events to grow, 
eventually causing total rupture of the fault jlBl . 
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FIG. 2. The frequency distribution for the number of top- 
plings and for the dimensionless event duration are shown for 
effective leakage factor <f> — 2.0; the plot with no leakage is 
shown for a comparison. 



This smoothing effect is reflected also in the distribu- 
tion for the number of topplings, shown in the lower plot 
on Fig. ||. Here, a counter is updated each time a site 
fails, even if it had failed before in the same avalanche. 
This number reflects more closely the model's equivalent 
for the moment release in an event. The fraction of mul- 
tiple failures vanishes in the exact mean-field limit, and 
the number of topplings become equivalent to the event 
area. This is no longer true for the model with stress 
leakage. The single maximum at the upper end of Fig. 
[j] corresponds, in fact, to very different moment releases, 
as shown in Fig. 0. An intriguing feature of the latter is 
the double bump at the large end, located close to inte- 
ger multiples of the system size. They suggest that total 
rupture of the system is more likely, in this region, than 
one would naively expect. 

The upper part of this plot shows the distribution of 
event durations, defined in the model as the number of 
updates in the fast time scale that are required for relax- 
ation. Again, the model with leakage reflects an excess 
of longer events over scaling: events with some range of 
large durations are more frequent, and the plot shows a 
local maximum close to its high end. 

The SCCA models with no leakage dynamics show no 
signs of a characteristic-event regime. The power spec- 
trum of the distribution of interevent times for large size 
events is white. The inclusion of leakage dynamics, how- 
ever, changes this aspect radically, as shown in Fig. [|. 
The distribution of interevent times has a pronounced 
maximum, corresponding to a characteristic period be- 
tween large rupture events. This feature is tuned by the 
leakage parameter, and will be more dominant as it in- 
creases. 
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FIG. 3. The frequency distribution for dimensionless in- 
terevent times between large events that rupture a fraction of 
0.9 of the total sites is shown for the dimensionless leakage 
factor <j} — 2.0. The plot shows the establishment of a charac- 
teristic event regime, superimposed on a complex background. 



The effects of the leakage dynamics on the statistical 
properties of the model's stress field are made more ex- 
plicit in Fig. |J. This plot shows the time evolution, in 
the slow time scale, of its average value and roughness, 
together with the configuration entropy. This last quan- 
tity is a measure of the degree of ordering of the stress 
distribution H] . Notice that the average stress is an in- 
creasing function of the time between large events, both 
for the model without and with leakage, but that in the 
latter the overall average is higher. As opposed to what 
was seen in a SCCA model with no leakage but with 
varying dissipation and weakening of failed sites |3 , we 
do not recognize a mode switching dynamics as present 
in our case. It remains to be seen what features would 
result from a combination of these dynamics. 

We presented in this paper a SCCA model with aseis- 
mic creep superimposed on a seismic threshold dynam- 
ics, as recent laboratory experiments have shown to ex- 
ist on solid-on-solid friction. In a mean-field approxi- 
mation, with infinite range of uniform interactions, the 
equations for the interevent stress time evolution can be 
solved. Computer simulations of the resulting SCCA 
model, combining the two dynamics, show a superpo- 
sition of complex time and space patterns with a more 
regular occurrence of characteristic events, with rupture 
of the entire fault. 
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FIG. 4. The plot compares the time evolution of the di- 
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